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Report  Title 

Normalized  Implicit  Radial  Models  for  Scattered  Point  Cloud  Data  without  Normal  Vectors 

ABSTRACT 

We  describe  some  new  methods  for  obtaining  a  mathematical  representation  of  a  surface  that  approximates  a  scattered 
point  cloud,  {(xyz)iN}  i  i  i , ,  =  1,L,  without  the  use  or  need  of  normal  vector  data.  The  fitting  surface  is  defined 
implicitly  as  the  level  set  of  a  field  function  which  is  a  linear  combination  of  trivariate  radial  basis  functions.  Optimal 
approximations  are  based  upon  normalized  least  squares  criteria  which  lead  to  eigenvalue/eigenvector 
characterizations.  The  normalized  aspect  allows  for  the  exclusion  of  the  need  of  normal  vector  estimates  which  is  one 
of  the  unique  features  of  this  new  method.  Localizing  techniques  are  introduced  to  allow  for  the  efficient  application  of 
these  new  methods  to  large  data  sets.  The  use  of  a  variety  of  radial  basis  functions  are  introduced  through  various 
examples  that  illustrate  the  performance  and  efficiency  of  the  new  methods. 


Normalized  Implicit  Radial  Models  for  Scattered  Point  Cloud  Data 

without  Normal  Vectors 


Gregory  M.  Nielson,  Tempe 
Abstract 

We  describe  some  new  methods  for  obtaining  a  mathematical  representation  of  a  surface 
that  approximates  a  scattered  point  cloud,  {(x; ,  y.  ,zr)  i  =  1 ,  •  •  ■ ,  /V }  without  the  use  or  need 

of  normal  vector  data.  The  fitting  surface  is  defined  implicitly  as  the  level  set  of  a  field 
function  which  is  a  linear  combination  of  trivariate  radial  basis  functions.  Optimal 
approximations  are  based  upon  normalized  least  squares  criteria  which  lead  to 
eigenvalue/eigenvector  characterizations.  The  normalized  aspect  allows  for  the 
exclusion  of  the  need  of  normal  vector  estimates  which  is  one  of  the  unique  features  of 
this  new  method.  Localizing  techniques  are  introduced  to  allow  for  the  efficient 
application  of  these  new  methods  to  large  data  sets.  The  use  of  a  variety  of  radial  basis 
functions  are  introduced  through  various  examples  that  illustrate  the  performance  and 
efficiency  of  the  new  methods 

Key  Words:  Surface  fitting,  point  clouds,  implicit  least  squares,  isosurfaces,  noisy  3D 
data,  scattered  data  approximation 


1.  Problem,  Introduction  and  Motivation 

We  present  a  new  technique  for  modeling  a  scattered  point  cloud  data, 
{(x, ,  v, ,  zt )  i  =  1  Based  upon  least  squares  error  criteria  the  method  determines 

the  parameters  of  an  implicit  model,  F(x,y,z),  so  that  the  zero  level  contour  surface, 
{(x,y,z):  F(x,y,z)=  0}  yields  a  surface  that  approximates  the  point  cloud.  One  of  the 
unique  features  of  this  new  method  is  the  lack  of  the  need  of  normal  vector  estimation. 
Estimating  normal  vectors  (including  orientation)  can  be  problematic  and  error  prone 
and  so  it  is  rather  beneficial  if  this  process  can  be  avoided.  Prior  to  proceeding  to  the 
details  of  describing  our  new  method,  we  give  a  brief  overview  of  point  cloud  fitting 
techniques  which  establishes  the  context  and  provides  motivation. 

One  of  the  first  things  that  we  should  point  out  is  that  the  problem  of  point  cloud  fitting 
should  be  distinguished  from  that  of  scattered  data  modeling  [11,  12,  21].  Even  though 
many  of  the  basic  techniques  and  tools  from  CAGD  (Computer  Aided  Geometric  Design) 
and  multivariate  approximation  theory  apply  to  both  problems,  they  are  basically 
different.  The  problem  of  scattered  data  modeling  is  concerned  with  methods  of 
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producing  a  bivariate  function  F(x,  y)  such  that  F(xi,yi)~  zt .  That  is,  for  the  traditional 

scattered  data  modeling  there  is  the  assumption  that  the  data  consists  of  samples  taken 
from  the  surface  graph  of  a  bivariate  function.  One  fundamental  connection  between  the 
two  problems  is  through  some  type  of  parameterization  of  the  point  cloud  [8,  9]  whether 
this  be  implicit  or  explicit.  The  tenn  “scattered  data”  was  coined  by  Schumaker  in  his 
1976  paper  and  there  was  a  great  deal  of  interest  and  published  research  on  (mainly) 
bivariate  problems  in  the  70s  and  80s.  With  the  advent  of  scientific  visualization  along 
with  volume  visualization  in  the  90s,  there  was  growing  interest  in  trivariate  scattered 
data  modeling  [21]  and  interest  in  this  area  continues  to  grow.  In  many  respects,  the 
problem  of  point  cloud  fitting  is  more  difficult  because  it  is  less  understood,  but  the 
widespread  and  strong  need  (see  [20])  for  practical  and  effective  methods  make  this  an 
important  problem. 

A  number  of  methods  involve  the  signed  distance  function  (see  [14]),  d(p),  which  is  a 
trivariate  function  defined  to  be  zero  on  the  surface  S ,  negative  interior  to  S  and 
positive  outside  of  S .  The  surface  is  extracted  from  d(p)  as  a  triangular  mesh  surface 
approximation  to  the  zero  level  isosurface.  Typically,  d(p)  is  sampled  on  a  3D 
rectilinear  grid  and  a  method  like  the  marching  cubes  algorithm  [22]  is  used.  Once  it  is 
decided  what  the  metric  or  the  definition  of  distance  from  a  surface  to  a  point  cloud  is  to 
be,  it  is  usually  not  too  difficult  to  develop  algorithms  for  the  efficient  computation  of  the 
distance  function.  The  particular  difficulty  here  is  getting  the  sign  right;  that  is,  to  be  able 
to  efficiently  and  effectively  detennine  when  a  point  is  inside  the  surface  or  outside. 
Typical  of  the  methods  based  upon  distance  functions  is  that  of  [15],  where  the  sign  is 
based  upon  local  least  squares  estimates  of  the  normal  vector  of  the  surface  and  a 
consistent  orientation  (in  or  out)  is  sought  with  the  Riemannian  map  estimate.  One  of  the 
drawbacks  to  this  method  is  the  heuristics  of  the  signed  distance  function  calculation  may 
lead  to  gaps  in  the  surface  and  the  difficulty  of  choosing  the  proper  resolution  for  the 
marching  cubes  voxel  grid  can  have  detrimental  effects  on  the  success  of  the  method. 

Another  important  concept  involved  in  point  cloud  fitting  problems  is  the  Delaunay 
tetrahedrization  and  its  dual,  the  Dirichlet  tessellation  and  Voronoi  diagram.  The  method 
based  upon  alpha  shapes  of  [6]  is  a  typical  and  early  example.  Here  the  first  step  is  the 
Delaunay  tetrahedrization.  The  second  step  is  to  apply  the  alpha-erasure  to  remove 
tetrahedra,  triangles  and  edges  whose  minimum  surrounding  sphere  is  not  contained  in 
the  alpha  -erasure  sphere.  The  result  is  called  the  alpha-shape.  In  the  third  step,  triangles 
for  the  final  surface  are  selected  so  that  a  sphere  of  radius  alpha  containing  the  triangle 
does  not  contain  any  other  point  cloud  points.  The  main  negative  aspect  of  this  approach 
is  the  choice  of  a  suitable  value  of  alpha.  Too  big  of  a  choice  leads  to  poor 
approximations  not  utilizing  many  of  the  points  of  the  point  cloud  and  too  small  a  choice 
leads  to  gaps  and  fragmented  surfaces.  We  mention  another  potential  drawback  to  these 
types  of  methods  for  certain  applications.  The  resulting  triangular  mesh  surface  has 
vertices  that  are  points  of  the  original  point  cloud.  For  noisy  data  or  overlapping  data 
resulting  from  imperfect  registration  of  scanned  data,  this  may  be  undesirable.  Rather 
than  interpolated  the  point  cloud  (or  a  subset),  it  is  potentially  more  desirable  to 
approximate  it  for  some  applications. 
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Another  class  of  methods  obtains  a  triangle  mesh  surface  (TMS)  approximation  to  the 
point  cloud  by  deforming  an  initial  approximation  using  metrics  that  measure  the  distance 
between  the  point  cloud  and  the  TMS  with  possible  refinement  of  the  TMS. 
Representative  examples  of  such  methods  include  [7]  and  [24].  Adaptive  techniques 
allow  fitting  methods  to  efficiently  apply  more  resources  to  regions  of  higher  complexity 
(see  [1],  [2]). 

2.0  Normalized  Implicit  Eigen-value/vector  Least  Squares  Models 

Given  a  collection  of  scattered  points  Pi,i  =  l,---,N  we  are  interested  to  find  a  field 
function  F(p)  of  an  implicitly  define  model  so  that  zero  level  contour, 
{(x,y,z):  F(x,y,z)  =  0},  represents  an  approximation  to  the  data  in  the  sense  that 
F(Pi ) «  0,  /  =  1,  •  •  • ,  A  .  We  assume  that  the  fitting  function  F  is  spanned  by  the  basis 
functions,  B  .(/*),  /  =  1 ,  ■  •  ■ ,  M  and  that  the  least  squares  criteria  is  used  and  so  we 
consider  the  problem: 


mm 
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Without  additional  constraints,  this  problem  is  not  well  formulated  since  an  optimal  fit  is 
the  identically  zero  function.  We  add  the  additional  constraint  of  nonnalizing  the 
weights/coefficients  of  the  fitting  function  and  consider  the  modified  optimization 
problem: 
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We  now  set  out  to  find  the  characterizing  equations  for  the  minimizing  set  of  coefficients, 
F.  /  =  1 ,  •  •  • ,  M  .  While  it  is  possible  to  invoke  some  general  results  concerning  Rayleigh 

quotients  (see  [16])  it  is  both  interesting  and  instructive  to  proceed  using  conventional 
calculus  techniques.  To  this  end,  we  introduce  the  following  notation  for  the  objective 
function  of  (1), 
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A  minimizer  of  (f> ,  necessarily  must  be  a  stationery  point;  that  is,  a  root  of  the  gradient, 
which  requires  that 


3(#(fi,fj,-,fJ'| 


v  ^(f,  ,  f3  ,  •  •  • ,  Fm  )  =  0  = 
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Using  the  standard  calculus  quotient  rule  for  derivatives,  we  have 
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j= i 

It  is  easy  to  see  that  the  above  system  of  equations  is  equivalent  to 
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where  A  =  BB*  denotes  the  gram  matrix  and 

' Bi(pi)  *ite)  -  bApnY 
B=  b2{px)  B2{P2 )  -  B2{Pn) 

kbm{px)  bm(p2)  ■■■  bm{pn\ 
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At  this  point,  we  assume  that  the  data,  Pi,i  =  l,--,N  and  the  basis  functions 
Bi,j  =  \,---,M  are  such  that  A  is  positive  definite.  From  Equation  (4),  we  can  see  that 

minimizing  solution,  F  must  be  an  eigenvector  of  A .  If  we  let  A  be  the  associated 
eigenvalue  then  we  may  also  conclude  that 
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and  so  we  have  established  the  fact  that  the  vector  F  that  minimizes  (2)  is  the 
eigenvector  of  the  gram  matrix,  BB*  associated  with  its  smallest  eigenvalue. 

2.1  Proof  of  Concept  Example 

We  now  present  a  very  simple  example  which  illustrates  the  fitting  process  and  which 
can  be  subsequently  used  by  others  to  check  and  verify  concepts  and  implementations. 
Here  we  have  thirteen  (13)  data  points:  (0.44,0.13),  (0. 24, 0.15),  (0.22,0.35), 
(0.31,0.49),  (0.35,0.60),  (0.43,0.70),  (0.51,0.77),  (0.57,0.70),  (0.70,0.55), 

(0.85,  0.30),  (0.88,  0.23),  (0.80,  0.10),  (0.60,  0.10)  which  are  depicted  in  the  right  image 
of  Figure  1.  We  use  a  set  of  basis  functions  which  give  rise  to  the  general  model 

F{x,  y)  =  a  +  bx  +  cy  +  ^  Fye  ^  Qj  ^ 

j= i 

where  P  =  {x,y)  and  =  are  the  “knots”  for  the  exponential  radial  basis 

functions.  For  this  particular  example,  we  simply  choose  R  =  1 ,  ju  =  l  M  =  3  and  knot 
values  Qx  =  (0.5 120,  0.6640) ,  Q2  =  (0.7825,  0.1825)  and  Q3  =  (0.3025,  0.2800).  These 

knots  are  selected  according  to  the  “venetia  criteria”  which  is  explained  below,  but 
basically  this  means  that  these  three  knots  are  the  “closest”  three  points  to  the  data  points 

{ ^  ]' = i  •  This  requires  that  each  knot  be  the  centroid  of  the  data  points  lying  in  its 
Voronoi/Dirchlet  cell/tile.  For  this  particular  case,  we  have 

Ql=(P5+P6  +  P1+Ps+Pg)/5,  Q2  =(Pl0  + Pn+ Pn+ Pu)/4,  Q3  =  (Px  +P2  +P3  +P4)/4 
and 
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The  smallest  eigenvalue  of  BB*  is  0.00273142  leading  to  the  solution 


F(x,y)  =  0.740  +  0.088x  -  0.093y  -  ,258eHl^a|  -  0.526eH|/"e21  -  0.306eH|i"a 
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Which  is  depicted  in  Figure  1. 
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Figure  1.  The  normalized  implicit  field  function  is  shown  on  the  left  and  on  the  right  is  shown  the  zero 
level  contour  of  this  implicit  fit  along  with  the  input  scattered  point  cloud. 

3.  Venetia  Criteria  for  Knot  Selection. 

Some  related  methods  for  knot  selection  have  previously  been  used  and  discussed  in  [10], 
[13],  [23],  [25].  This  method  of  selecting  knots  used  here  is  based  upon  minimizing  the 
distance  between  the  knots  {Q  ,j  =  and  the  scattered  point  cloud, 

{  P;,i  =  This  minimal  distance  requirement  leads  to  the  sufficient  condition 

called  the  “venetia  criterion”  where  Qf  is  required  to  be  the  centroid  of  the  points  lying 

in  its  Theissen/Dirichlet  region.  The  Theissen/Dirichlet  region  is  defined  to  be  the  set  of 
points  closer  to  Qjt  han  any  other  Qk ,  k  ^  j  and  is  denoted  by  T /£>[0;J.  Computational 

algorithms  for  the  case  of  two-dimensions  are  discussed  in  [10]  and  for  case  of  three- 
dimensions  in  [25].  These  algorithms  involve  the  basic  “venetia  iteration”  where  a  knot 
QP  is  updated  to  a  new  knot  (T  AVI  i  which  is  the  centroid  of  the  data  points  lying  in  the 

Theissen/Dirichlet  region  of  QP  .  As  it  is  pointed  out  in  both  [10]  and  [25],  this  type  of 

iteration  quickly  converges,  but  not  necessarily  always  to  a  global  optimal  knot  set  that  is 
a  minimal  distance  to  the  data  points.  Nevertheless,  we  have  found  that  with  good  initial 
knot  distribution,  this  type  of  iteration  leads  to  a  very  efficient  means  of  determining  knot 
locations  which  subsequently  leads  to  overall  efficient  implicit  approximation  models. 

In  the  example  shown  in  Figure  2,  we  use  as  an  initial  approximation  the  centroids  of  the 
points  lying  in  the  lattice  grid  points  of  a  3x3  grid.  As  it  turns  out  two  of  the  nine  cells 
are  void  of  data  points  and  so  this  approach  only  leads  to  seven  knots. 
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Figure  2.  The  scattered  point  cloud  is  illustrated  with  “unfilled”  circles  and  the  knots  are  the  “filled” 
circles.  The  knots  in  the  right  image  are  characterized  by  the  “Venetia  Criterion”  which  is  a  sufficient 
condition  for  minimizing  the  overall  distance  between  the  knots  and  the  scattered  point  cloud  and  implies 
that  each  knot  is  the  centroid  of  the  data  points  lying  in  their  respective  Theissen/Dirichlet  regions.  Each  of 
the  nine  regions  of  the  left  image  containing  data  points  is  “seeded”  with  a  knot  prior  to  the  optimization 
process. 


4.0  The  Localization  Technique 


For  large  data  sets,  the  concepts  of  localization  can  be  used  where  the  fitting  process  is 
broken  down  into  s  a  collection  of  smaller  problems  and  the  final  model  is  a  blend  of 
these  local  smaller  and  simpler  fitting  problems.  This  general  technique  is  described  in 
[21].  The  basic  ideas  of  the  technique  were  also  used  by  [26]  where  they  refer  to  Franke 
&  Nielson  [11]  as  their  original  source  for  the  basic  idea.  Here,  we  give  a  brief,  high 
level,  overview  of  this  technique  within  the  context  of  scattered  data  interpolation.  Input 
to  the  problem  of  scattered  data  modeling  consists  of  a  collection  of  scattered  data  point 
locations  Pt,i  =  !,•••, TV  contained  in  the  domain  D  and  associated  dependent  data  values 


Yt,i  =  .  It  is  desired  to  construct  an  interpolating  function  F  defined  over  D  that 

has  the  property,  F(Pj)=Yj,i  =  1  ,--N.  Let  Wk,k  =  be  a  collection  of  functions 

defined  overD  and  further  assume  that  the  each  Wk  has  local  support.  That  is  the  set  of 
points  Pj  such  that  Wk  (Pf )  ^  0  is  a  small  subset  of  the  total  number  of  N  points  and  the 
every  data  point  must  lie  in  the  support  of  some  Wm .  Also  it  must  be  the  case  that 

M 

Ywk(r)  =  1  for  all  points  P  in  the  composite  domain.  Note  that  this  last  condition  is 

k=l 


not  a  deal  breaker  for  if  it  is  not  satisfied,  then  it  is  possible  to  use  the  “normalized’' 


weight  functions  wk  = 


W, 


IX 


which  do  have  this  sum-to-unity  property.  Assume  we 


have  a  collection  of  local  interpolating  functions  Fk(P )  each  having  the  property  that 
fAp,)=Y,  for  all  P,  in  the  support  of  Fk(p),  then  the  composite  weighted  model 


M 

F(p)  =  ^  Wk{P)Fk(P )  has  the  property  that 

k=\ 


7 


F  (P, ) =  S  w>  (p,  K  (p) =  r<  •  for  a11  >'  =  1. ' '  ■  • .  ■ ff  ■ 

k= 1 


Here  we  are  not  solving  an  interpolation  problem,  but  rather  a  least  squares 
approximation  and  also  we  are  not  doing  conventional  scattered  data  interpolation  (  i.  e. 
Fk{p)=Yi  )>  but  yet,  we  can  still  exploit  these  basic  ideas  to  make  large  data  set 
problems  tractable  with  our  present  approach.  We  can  explain  best  how  this  technique 
works  with  a  simple  example  with  32  scattered  data  points  as  illustrated  in  Figure  3.  We 
wish  to  find  an  implicit  model  similar  to  what  was  accomplished  in  the  previous  example 
of  Figure  1,  but  rather  than  solve  the  overall  problem,  we  break  it  down  into  some  smaller 
problems  and  combine  the  solutions  using  the  techniques  of  localization.  Potential 
candidates  for  localizing  functions  are  piecewise  bilinear/trilinear  functions  defined  over 
a  user  specified  NxM  or  NxM  x  K  grid.  Each  localizing  function  Wtj  has  the 

property  that  wJk,l)  =  8ik8 7 .  Since  the  sum  of  these  localizing  functions  is  a  piecewise 

bilinear  (trilinear)  which  takes  on  the  value  1  at  each  grid  point  and  the  summation  is 
identically  equal  to  1.  For  each  localizing  function,  we  compute  the  approximation  based 
upon  the  data  points  lying  in  the  support  of  this  localizing  function.  The  final 
approximation  is  then  the  weighted  sum  of  the  local  approximations.  For  the  example  of 
Figure  3,  the  domain  is  the  square  [0,3]x[0,3].  Actually  for  this  small  case  where  the 
“exceptional”  boundary  cases  constitute  such  a  large  percentage  of  the  cases,  we  only 
compute  four  localized  approximations;  one  each  of  the  four  subsets  of  data  lying  in  the 
domains  [0,2]x[0,2],  [0,2]x[l,3],  [l,3]x[0,2]  and  [l,3]x[l,3].  These  are  respectively 
referred  to  as  Fn,F12,F2l  and  F22 .  The  final  result  is 

G(x,  y )  =  (w00  +  w01  +  w10  +  w, ,  )Ft ,  (x,  y)  +  (w02  +  w03  +  w12  +  w13  )f12  (x,  y ) 

(w20  +  W21  +  w30  +  w31  )F2l  (x,  y)  +  (w22  +  w23  +  w31  +  w33  )F22  (x,  y ) 

which  is  illustrated  in  the  two  lower  images  of  Figure  3.  Using  only  four  approximations 
here  rather  than  16  serves  to  point  out  the  possibility  of  other  data  dependant 
consolidations  to  be  used  for  localization  that  are  based  upon  rectilinear  grids.  There  are 
a  large  number  of  potentially  interesting  and  useful  options  to  consider  and  exploit  in  this 
context. 


Figure  3.  An  example  that  illustrates  the  localization  process.  The  upper  four  graphs  show  the  data  and  the 
contour  of  the  normalized  implicit  model  applied  to  the  data  contained  in  bounded  rectangular  region 
respectively.  These  regions  constitute  the  support  of  the  localizing  functions  being  used.  The  two  images 
in  the  bottom  row  show  the  final  blend  of  these  four  local  approximations.  The  blending  functions  are 
piecewise  bilinear  functions  which  take  on  0  or  1  at  each  of  the  integer  lattice  domain  points  and  have 
support  on  the  respective  rectangles.  The  lower,  left  image  is  the  global  normalized  implicit  model  and  the 
lower,  right  image  shows  the  zero  level  contour  along  with  the  input  scattered  point  cloud.  It  is  really  very 
interesting  that  the  quality  of  the  fit  is  so  good  even  without  the  use  of  normal  vector  (orientation) 
information. 


5.  Examples  and  Applications 

5.1  Parameter  specification  for  least  squares  parametric  curve  fitting 

This  first  example  is  a  typical  representative  of  a  2D  curve  fitting  problem  where  we  are 
given  a  collection  of  points  (xi,yl),i  =  \,...,N  and  it  is  desired  to  fit  a  curve  that 
approximates  the  shape  inferred  by  these  points.  Since  the  data  does  not  necessarily  infer 
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a  function  relationship,  y  =  F(x),  parametric  methods  (x(t),y(t ))  0  <  t  <  1  can  be  used. 
Most  any  approach  based  upon  parametric  curves  will  require  the  specification  of 
associated  parameter  values  ti  for  each  scattered  point  (x, ,  y. ) .  For  this  example  we  will 

use  least  squares  fitting  and  so  once  the  parameter  values  are  specified,  the  linear,  least 
squares  fitting  problem  is  rather  straight  forward: 


N  f  M 
ic*'  i=i\k=i 


(V) 


where  {BvB2,--,Bm}  are  the  basis  functions,  p;  =(xj,yi)  are  the  points  of  the  scattered 
point  cloud,  ck,k  =  are  the  to  be  detennined  coefficients  of  the  final  fitting  model 

and  t.,i  =  are  the  all-important  parameter  values  which  have  yet  to  be  specified. 

Once  these  parameter  values  are  specified,  solving  the  minimization  problem  of  (7)  has 
been  studied  extensively  and  there  are  many  efficient  and  effective  computational 
methods  generally  available.  If  these  parameter  values  are  not  known,  but  the  scattered 
data  points  are  at  least  ordered,  then  there  are  available  a  number  of  rather  effective 
means  for  determining  parameter  values.  If  we  choose  to  not  specify  these  values,  then, 
it  is  possible  that  these  values  can  be  left  undetennined  and  simply  be  added  to  the 
unknowns  for  the  minimization  process.  This  makes  the  least  squares  problem  nonlinear 
and  extremely  problematic.  Here,  we  suggest  using  the  implicit  methods  as  a  basis  for 
selecting  these  parameter  values.  In  a  nutshell,  the  approach  is  as  follows.  First  an 
implicit  model  as  we  have  described  here  is  fit  to  the  data  points  (x; ,  y. ) .  Next  the  data 

points  are  projected  onto  the  contour  of  the  implicit  model.  These  values  are  ordered  by 
there  projections  onto  the  zero  level  contour.  This  ordering  (and  possibly  the  relative 
distances  along  the  contour)  is  used  to  detennine  the  parameter  values.  In  the  example  of 
Figure  4,  we  use  an  approximate  projection  of  each  data  point  onto  the  contour  of  the 
implicit  model.  The  projection  is  taken  to  be  the  intersection  of  the  contour  and  the 
parametric  line 


L{t)  =  {xiyl) 


+  t 


dF  (  ^dF 
^{xiyi),— 
ox  oy 


Depending  upon  the  details  of  how  F  is  stored  and/or  represented,  the  gradient  may  or 
may  not  be  immediately  available.  If  not,  discrete  approximations  can  be  used. 

For  the  example  of  Figure  4,  we  have  70  scattered  data  points.  An  implicit  model, 
similar  to  that  of  Figure  3  is  fit  to  these  points.  The  data  and  the  fit  are  shown  in  the  top 
image  of  Figure  4.  Each  of  the  data  points  is  projected  onto  the  contour  of  the  implicit 
model.  These  projected  points  give  a  relative  order  for  the  data  points  which  is  used  to 
obtain  the  associated  parameter  values  }™  .  Based  upon  these  parameter  values,  a 

periodic,  parametric  cubic  spline  is  fit  to  the  airfoil  data  in  the  least  squares  sense.  The 
result  is  shown  in  the  bottom  image  of  Figure  4. 
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Figure  4.  Top  image  shows  airfoil  data  and  NIRM  fit.  The  bottom  image  shows  the  least  squares 
parametric  spline  fit  with  proportionately  spaced  knots  so  that  each  knot  interval  has  approximately  the 
same  number  of  data  point  parameter  values. 

5.2  Boot  and  Foot  CSG  Example 

This  next  example  utilizes  a  scattered  point  cloud  data  set  obtained  from  scanning  a 
physical  boot  which  is  depicted  in  the  left  image  of  Figure  5.  There  are  approximately 
one  million  data  points.  While  the  point  set  is  fairly  uniformly  dense,  it  is  rather  noisy 
due  to  the  method  of  collection  which  requires  only  a  video  projector  and  a  camera  (see 
[27]  for  example).  A  10%  sample  of  the  points  is  shown  in  the  center  image  of  Figure  5. 
A  rendering  of  the  implicit  model  is  shown  in  the  far  right  image  of  Figure  5.  A  localized 
model  using  the  techniques  of  Section  4  is  used.  The  grid  analogous  to  the  3x3  grid  of 
Figure  2  is  a  uniform  grid  of  size  100x100x50  over  the  bounding  cuboid  of  the  data.  The 
localizing  technique  of  Section  4  is  used  based  upon  piecewise,  trilinear  localizing 
functions  defined  over  a  50x50x25  grid.  The  basis  functions  are  of  the  form 

M 

f(x,  y,  z)  =  a  +  bx  +  cy  +  dz  +  If 

j= i 

These  basis  functions  have  previously  been  successfully  used  in  [23].  The  value  of  R  is 
the  radius  of  the  support  regions  of  the  localizing  functions.  The  power  /u  is  simply  set 
to  2  for  this  particular  example.  The  knot  selection  technique  described  above  is  seeded 
with  7  knots  per  non-empty  cell.  The  final  implicit  fit  is  sampled  on  a  3D  rectilinear  grid 
and  the  dual  marching  cubes  [22]  method  is  used  to  compute  a  polygon  contour  surface. 


rj-F-Qj 

R 


f'j 
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Figure  5.  Left  image  is  photograph  of  motocross  boot.  Middle  image  is  rendering  of  10%  of  the 
approximately  one  million  data  points  scanned  from  a  physical  boot.  Right  image  is  a  rendering  of 
isosurface  extracted  from  implicit  model. 


Figure  6.  Top  left  image  is  contour  of  implicit  model  of  “noisy  foot”  data  using  exponential  basis 
functions.  Bottom  left  is  right  foot  obtained  as  F(y,X,z)  =  O).  Right  image  illustrated  the  results  of 
Boolean  operations  which  are  easily  accomplished  with  implicit  models. 

In  Figure  6  we  show  yet  another  example.  This  example  utilizes  the  “noisy  foot”  data 
from  Geomagic.  Here  we  use  the  3D  versions  of  the  same  exponential  basis  functions 
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used  in  Example  2.1.  In  the  right  image  we  show  an  example  of  where  the  “foot”  has 
been  removed  from  the  “boot”.  One  of  the  reasons  for  including  this  example  is  to 
underscore  the  ease  of  computing  Boolean  operations  (see  [3])  on  models  that  are 
implicitly  defined  as  is  the  case  for  the  fitting  models  of  this  paper.  If  we  let  F A  denote 
the  field  function  for  a  point  set  A  whose  boundary  is  a  surface  of  interest.  That  is 
A  =  {(x  ,y,z ) :  F A  {x,  y,z)>  0} .  And  if  B  is  another  three  dimensional  point  set  with  the 
field  function  F B ,  then  it  is  easy  to  see  that  the  union  is  defined  by 

A  u  B  =  { {x,  y,  z) :  Max{FA  (x,  y,  z),  FB  (x,  y,  z))  }  (8) 

and  so  we  have  that  FAjH  =  Max{F A,FB).  Similarly,  we  have  for  the  intersection 

FAnB  =Mw{Fa,Fb). 


6.  Remarks 


1 .  Here  we  have  used  two  classes  of  radial  basis  functions.  Namely 


and 


M 


F{x  ,y,z)  =  a  +  bx  +  cy  +  dz  +  IV 

j= i 


-R,\p-Qi 


(9) 


F[x,y,z)  =  a  +  bx  +  cy  +  dz  +  ^  F  f 


Rj-P-Qj 

R: 


(10) 


There  are  a  great  number  of  other  radial  basis  functions  which  could  be  considered  that 
may  have  particularly  efficient  fitting  capabilities  for  certain  types  of  data  sets. 


2.  Using  the  same  data  as  in  the  “proof  of  concept”  example  of  Section  2.1,  but  with  the 
basis  functions  of  (TO)  (with  R ,  =  1 ,  =  2  v  j),  we  obtain  the  results  shown  in  Figure 

7.  While  the  contour  seems  to  be  a  superior  fit  to  that  of  the  earlier  example,  the  criteria 
of  normalized  least  squares  is  actually  larger  for  this  model.  The  smallest  eigenvalue  for 
example  of  Figure  1  is  .00273  and  here  the  smallest  eigenvalue  is  0.00767  and  so,  in  a 
standard  least  squares  sense  of  measuring  error  the  approximation  of  Figure  1  is  actually 
better.  This  seems  to  indicate  that  some  possible  adjustment  in  the  overall  fitting  criteria 
may  be  in  order.  A  shallow  field  function  would  lead  to  a  small  RMS  error,  but  the 
distance  of  the  point  cloud  to  the  contour  surface  could  still  not  be  small.  Normalizing 
with  gradient  norms  is  a  possibility,  but,  of  course,  we  don’t  want  to  destroy  the 
relationship  of  the  solution  we  have  here  and  its  connection  to  the  eigenvalues  and 
eigenvectors.  Care  must  be  taken  to  appropriately  alter  the  fitting  norm,  but  still  be  able 
to  solve  the  problem  without  resorting  to  the  iterative  methods  of  nonlinear  least  squares 
fitting.  We  will  report  on  our  progress  in  this  area  in  a  future  paper. 


13 


Figure  7.  The  same  point  cloud  as  in  the  example  of  Figure  1,  but  with  different  radial  basis  functions. 
The  contour  appears  to  be  a  better  approximation,  yet  the  implicit  least  squares  fit  error  is  larger.  This 
motivates  discussion  about  using  a  different  fitting  criteria  (error  metric).  The  triad  of  the  right  image  is  the 
Dirichlet  tesselation  of  the  knots  which  are  indicated  with  numeric  values.  We  note  that  the  knots  satisfy 
the  “ventia”  criteria  in  that  they  are  the  centroids  of  the  scattered  point  cloud  lying  in  the  respective 
Dirichlet  tiles. 
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